Performance of High Order Filter Methods for a Richtmyer-Meshkov Instability

Sixth-order compact and non-compact filter schemes that were designed for multiscale Navier-Stokes, and ideal and non-ideal magnetohydrodynamics (MHD) systems are employed to simulate a 2-D Rightmyer-Meshkov instability (RMI). The suppression of this RMI in the presence of a magnetic field was investigated by Samtaney (2003) and Wheatley et al. (2005). Numerical results illustrated here exhibit behavior similar to the work of Samtaney. Due to the different amounts and different types of numerical dissipations contained in each scheme, the structures and the growth of eddies for this chaotic-like inviscid gas dynamics RMI case are highly grid size and scheme dependent, even with many levels of refinement. 1 Numerical Method and Objective Methods commonly used for shock/turbulence interactions relying on switching between spectral or high order compact schemes and shock-capturing schemes are not practical for multiscale shock interactions. Frequent switching between these two types of schemes can create severe numerical instability. Our highly parallelizable class of high order filter schemes [7, 5, 8, 9, 10, 12, 13] does not rely on switching between schemes to avoid the related numerical instability. They have built-in flow sensors to control the amount and types of numerical dissipation only where needed, leaving the rest of the flow region free of numerical dissipation. Instead of solely relying on very high order high-resolution shock-capturing methods for accuracy, the filter schemes take advantage of the effectiveness of the nonlinear dissipation contained in good shock-capturing schemes and standard linear filters (and/or high order linear dissipation) as stabilizing mechanisms at locations where needed. The method consists of two steps, a full time step of spatially high order non-dissipative 2 B. Sjögreen and H.C.Yee (or very low dissipative) base scheme and an adaptive multistep filter consisting of the products of wavelet based flow sensors and linear and nonlinear numerical dissipations to filter the solution. The adaptive numerical dissipation control idea is very general and can be used in conjunction with spectral, spectral element, finite element, discontinuous Galerkin, finite volume and finite difference spatial base schemes. The type of shock-capturing scheme used as nonlinear dissipation is very general and can be any dissipative portion of a high resolution TVD, MUSCL, ENO, or WENO shock-capturing method. By design, the flow sensors, spatial base schemes and linear and nonlinear dissipation models are stand alone modules, and a whole class of low dissipative high order schemes can be derived at ease. The objective of this work is to illustrate the performance of our sixth-order low dissipation filter schemes for a 2-D inviscid Richtmyer-Meshkov instability (RMI) problem. 2 RMI Test Problem and Numerical Results RMI occurs when an incident shock accelerates an interface between two fluids of different densities. This interfacial instability was theoretically predicted by Richtmyer [3] and experimentally observed by Meshkov [2]. For the present study, the RMI problem investigated by Samtaney [4] and Wheatley et al. [6] as indicated in Fig. 1 has been chosen. The mathematical models are the 2-D Euler gas dynamics equations and the ideal MHD equations. The computational domain is −2 < x < 6 and 0 < y < 1. A planar shock at x = −0.2 is moving (left to right) toward the density interface with an incline angle of θ with the lower end initialized at x = 0. The density ratio across the interface is denoted by η, and the nondimensional strength of the magnetic field β = 2p0/B 0 , where the initial pressure in the preshocked regions is p0 = 1, and B0 is the initial magnetic field. The initial magnetic field is uniform in the (x, y) plane and perpendicular to the incident shock front. Numerical results shown below are forM = 2, θ = 45, η = 3, β−1 = 0 (Euler gas dynamics) and β−1 = 0.5 (magnetic field present). The computation stops at an evolution time t = 3.33. For this set of parameters and all studied numerical schemes, instability occurs near t = 1.8 for the gas dynamics case but not for the MHD case for the entire time evolution. Our numerical results exhibit behavior similar to the study of Samtaney. Computations by the sixth-order centered spatial compact base scheme with the compact linear filter [1], in conjunction with a second step nonlinear WENO5 filter (WENOfi) denoted by Comp66+Compfi+WENOfi using a 801 × 101 grid is shown in Fig. 2 (left) for the inviscid gas dynamics (GD) and the ideal MHD equations. Here WENOfi means the dissipative portion of the fifth-order WENO scheme (WENO5) in conjunction with our wavelet flow sensor as the nonlinear filter [10, 12, 13]. The same computation using the sixth-order central spatial (non-compact) base scheme in conjunction with High Order Filter Methods for RM Instability 3 WENOfi denoted by CEN66+WENOfi is shown in Fig. 2 (right). The classical fourth-order Runge-Kutta method is used for the sixth-order compact and non-compact filter schemes. For this low resolution grid, the accuracy between the two filter methods is similar. See Fig. 3 for the grid refinement study below. Computations using Comp66+WENOfi (i.e., without the linear compact filter step) or Comp66+Compfii (i.e., without the nonlinear WENOfi filter step) indicate spurious oscillations around shock regions. The present study arrives at the same conclusion drawn in [7, 12] on the behavior of compact spatial schemes for problems containing multiscale shock interaction. High order compact schemes are methods of choice for many incompressible and low speed turbulent/acoustic flows due to their advantage of requiring very low number of grid points per wavelength. In the presence of multiscale shock interactions, however, this desired property of high order compact base schemes seems to have diminished in both the gas dynamic and MHD test cases that we have studied (compared with the same order of accuracy of non-compact central base schemes). Also the compact spatial base scheme requires more CPU time per time step and it is less compatible with parallel computations than the central spatial base scheme. Consequently, the compact spatial base scheme requires added CPU time in a parallel computer framework. Figure 3 shows the inviscid gas dynamics comparison between CEN66+WENOfi and a second-order MUSCL for four grids (801× 101, 1601× 201, 3201× 401, 6401× 801). Here, computations using a second-order MUSCL and a secondorder Runge-Kutta method are denoted by MUSCL. Not shown are the same computation using CEN66 as the base scheme in conjunction with the dissipative portion of the Harten-Yee scheme and the wavelet flow sensor as the nonlinear filter (CEN66+HYfi). For similar resolution, MUSCL requires nearly 3 times finer grid size per spatial direction than CEN66+WENOfi and CEN66+HYfi. The eddy structures are different among the three filter methods and they are very different from the Samtaney adaptive mesh refinement (AMR) simulation with an equivalent uniform grid of 16384×2048. Note that except for WENOfi, the van Leer version of the van Albada limiter is used. For the second-order MUSCL scheme, the limiter is applied to the primitive variables. All methods employed the Roe’s approximate Riemann solver for the gas dynamics case and the Gallice approximate Riemann solver for the MHD case using our method of solving the conservative MHD system [9]. The present study can serve as a good example of failure of grid refinement for unsteady chaotic-like inviscid flow. As the grid is refined (in conjunction with different amounts and types of numerical dissipations contained in each scheme), smaller and smaller eddies are formed which combine to affect global flow through the energy cascade effect. Future work, including viscous effects and extensive comparisons is in progress [14]. 4 B. Sjögreen and H.C.Yee 3 Concluding Remarks The efficiency, accuracy and flexibility of the present class of low dissipative high order filter schemes rest on the fact that the multistep filter can be applied after a full time step of the Runge-Kutta method, and that a wavelet flow sensor is employed to control the amount of numerical dissipation where needed. The major CPU time intensive part of the shock-capturing computation is the nonlinear filter. In fact, the filter procedure requires slightly more CPU time than the Harten-Yee and MUSCL schemes. This is due to the fact that all sixth-order (or any order) filter schemes require only one Riemann solve per time step per direction (independent of the time discretizations of the base scheme step) as opposed to two Riemann solves per time step per direction by the MUSCL and Harten-Yee schemes using a second-order RungeKutta method. WENO5-RK4 requires at least twice the CPU time of all other methods since four Riemann solves per time step per direction are required by WENO5-RK4. RK4 stands for the classical 4th-order Runge-Kutta temporal discretization and WENO5 is the 5th-order WENO spatial scheme. Fig. 1. Problem definition High Order Filter Methods for RM Instability 5 Fig. 2. Comparison between Euler gas dynamics and MHD for the 6th-order compact spatial base scheme (left) and the 6th-order central (non-compact) spatial base scheme (right) using a (801×101) grid at t = 3.33. MHD solutions shown are mirror images of the original computations.